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Abstract 

In this paper, we consider signals with intra-wave frequency modulation. To handle this 
kind of signals effectively, we generalize our data-driven time-frequency analysis by using a 
shape function to describe the intra-wave frequency modulation. The idea of using a shape 
function in time-frequency analysis was first proposed by Wu in [R], A shape function could 
be any 27r-periodic function. Based on this model, we propose to solve an optimization problem 
to extract the shape function. By exploring the fact that s is a periodic function of 9 , we can 
identify certain low rank structure of the signal. This structure enables us to extract the shape 
function from the signal. To test the robustness of our method, we apply our method on several 
synthetic and real signals. The results are very encouraging. 


1 Introduction 

Nowadays, data play a more and more important role in our life. At the same time, the scientific 
community face a challenging problem: how to efficiently extract useful information from massive 
amount of data. In many real world problems, especially in engineering and physical problems, 
frequencies of the signal are usually very useful to help us understand the underlying physical 
mechanism. Hence, many time frequency analysis methods have been developed, for instance, the 
windowed Fourier transform, the wavelet transform mm, the Wigner-Ville distribution [8j, etc. 
In recent years, an adaptive time frequency analysis method, the Empirical Mode Decomposition 
(EMD) method [H [16] was developed. This method provides an efficient adaptive method to 
extract frequency information. The EMD method has found to be very useful in many applications. 
EMD method is purely empirical, it still lacks a rigorous mathematical foundation. Recently, 
several methods have been proposed attempting to provide a mathematical foundation for EMD 
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method, see e.g. the synchrosqueezed wavelet transform 0, the Empirical wavelet transform [B], 
the variational mode decomposition [9]. 

In the last few years, inspired by the EMD method and compressive sensing [1100, we proposed 
a novel time-frequency analysis method based on the sparsest time-frequency representation of 
multiscale data [TZi . In this method, the signal is decomposed into several components 

M 

m = £ a,j(t) cos Oj(t) + r(t), t € R, (1) 

3 = 1 

where aj(f), 9j(t) are smooth functions, 9'-{t) >0, j = 1, • • • , M, and r(t) is a small residual. We 
assume that a,j ( t ) and O'- are less oscillatory than cos 9j ( t ). The exact meaning of less oscillatory 
will be made clear later. We call aj(t) cos 9j(t) as the Intrinsic Mode Functions (IMFs) [J]. 

One main difficulty in computing the decomposition 0) is that the decomposition is not unique. 
To pick up the ’’best” decomposition, we proposed to decompose the signal by looking for the 
sparsest decomposition by solving a nonlinear optimization problem: 

Minimize M 

(ak)l<k<M ,(9k)l<k<M 

M 

Subject to: / = a*, cos 9k, 

k =i 

ak cos 9k € T>, 

where T> is the dictionary consist of all IMFs (see [12] for its precise definition). 

To solve ([2]), we proposed two algorithms. The first one is based on matching pursuit [123 
and the other one is based basis pursuit m ■ In a subsequent paper nn. the authors proved the 
convergence of their nonlinear matching pursuit algorithm for periodic data that satisfy certain 
scale separation property. 

Although the model (0) has been applied to a number of applications with success and the 
decomposition methods have been shown to be effective and efficient, there are also some applica¬ 
tions such as the Stokes waves or some nonlinear dynamic systems for which our methods do not 
work well. An essential difficulty for this type of data is that the instantaneous frequency, 9' k , is as 
oscillatory as or even more oscillatory than cos 9k, which we call intra-wave frequency modulation. 
One consequence is that the data may not have sparse representation in model (0). To effectively 
decompose the signal with intra-wave frequency modulation, we need another model to retrieve the 
sparse stucture. Our approach is to introduce a periodic shape function m to replace the cosine 
function in 0), then we get following model: 

M 

f{t ) = a k(t)s k {0k{t )) + r(t), (2) 

k=i 
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where s k is an unknown 27r-periodic ‘shape function’ and is adapted to the signal. The envelope 
a k (t) and the phase function 0 k (t) are smooth functions and are less oscillatory than cos 0 k (t) . We 
also assume that 0' k (t) > 0. 

With the introduction of shape function, even the signal with intra-wave frequency modulation 
has a sparse representation. But it is more difficult to find this sparse decomposition since the 
dictionary is larger than that in Jl]) due to the extra degree of freedom of shape functions. 

In this paper, we will introduce a method to extract the shape function for the signals with 
only one dominated shape function. First, the phase function is computed by our data driven 
time frequency analysis |12j . Once the phase function is obtained, by exploring the fact that s is 
a periodic function of 0 , we can identify certain low rank structure of the signal. This structure 
enables us to extract the shape function from the signal. 

The rest of this paper is organized as follows. In Section 2, the decomposition model for data 
with intra-wave frequency modulation is presented. The details of the algorithm and the localized 
version are given in Section 3 and 4. We present some numerical results in Section 5. Some 
concluding remarks are made in Section 6. 

2 Models for signal with intra-wave frequency modulation 

In order to design a computational algorithm for the model (|2|) . we first need to define the meaning 
of ’’less oscillatory”. With a given phase function 6(t), we construct a linear space V(0) which is 
spanned by the following basis 



where A < 1/2 is a parameter to control the smoothness of functions in V(0, A) and Lg = (0(1) — 
0(O))/27t is a positive integer and [0,1] is the span of the signal in time. And we require a k {t) and 
0 k (t) to be V(0 k , A) to enforce that they are less oscillatory than cos 0 k . 

Then, the model of the signal is given as follows: 

M 

f(t ) = ^ ~2a k (t)s k (0 k {t )), a k ,0' k &V(0 k ) 

k =1 

and 0' k > 0, s k is 27r-periodic. (4) 

Corresponding to this model, we propose to solve the following optimization problem to find 
the sparsest decomposition (pf|), 

min M (5) 

a k 

M 

Subject to: E ^k^k^fik) /) ®fc'Sfc(0fc) £ -Ad. 

fc=l 
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where the dictionary M is defined as 


M 




a k ,6' k eV(e k ), 9\. > 0 : 
s k is 27r-period function. 


( 6 ) 


This optimization problem is very difficult to solve. In this paper, we focus on a simpler case. We 
assume the signal is dominated by one component in M, i.e. 


m 


a(t)s(9(t)) + r(t), a,6'eV(9), 

and O' > 0, s is 27r-periodic. 


(7) 


Here, r(t) is the residual. The residual r(t) could be noise or trend or some kind of perturbation. 
No matter what r(t) is, we assume that it is small in amplitude compared with a(t)s(9(t)). 

Using the idea of matching pursuit, the decomposition in ([7]) can be obtained by solving the 
following optimization problem: 

min ||/(t)-a(i)s(0(t))||jj, (8) 

a,v,s 

subject to: a(t)s ( 6(t )) € A4. 

Although this optimization problem is much simpler than (JSJ), it is still very difficult to solve. It is 
highly nonlinear. The envelope a, the phase function 0, and the shape function s are all unknown. 
They are all adaptive to the data. This feature makes our method fully adaptive to the signal, but 
it also introduces additional difficulty to solve the resulting optimization problem ([8]) . 

Inspired by our previous work in the data-driven time-frequency analysis [T2|, we develop an 
efficient method to solve (|8|). First, the phase function is computed by our data-driven time- 
frequency analysis [12J. Once the phase function is obtained, by exploring the fact that s is a 
periodic function of 9, we can identify certain low rank structure of the signal. This structure 
enables us to extract the shape function from the signal. The details will be given in the next 
section. 


3 An efficient algorithm to compute the shape function 

First, to simplify the discussion, we assume that the phase function 9 has been obtained, then the 
optimization problem <[8|) can be reduced to the following problem: 

min || f(t)-a(t)s(0(t)) \\%, 

a,s 

subject to: a€V(9), s(-) is 27r-periodic. 

Since s is periodic, it can be represented by Fourier basis, 

K 

s(0) = £ c k e ike . (9) 

k=—K 
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In the above Fourier representation, we assume that s is iT-band limited, which is a good approxi¬ 
mation as long as s is smooth enough and K is large enough. 

Using the above representation, the optimization problem becomes 


min 

a,Ck 


K 

||/ —a ^ III, subject to: a € V{9). 

k=—K 


Next, in order to further simplify above optimization problem, we replace the standard l 2 norm in 
the objective function to l 2 norm in 0 space, 


2,0 — 



1/2 


where 8 = 6/Lg is normalized phase function which is used as the coordinate function in the 0-space 
and Lg = j n this paper, we assume that the signal lies in [0,1]. 

Then, the above optimization problem is reduced to 


+oo K z 

£ le^ e ' Le - a J2 c k e ik9 , 

u=—oo k=—K 2,0 

subject to: a£V(d). 

where fg is the Fourier coefficients of / in the 0-space, 

f g (u)= f 1 f{t)e-^e\t)dt. 

Jo 

Next, we represent the envelope a by the Fourier basis in the 0-space, 


min 

a,Ck 


Then we have 


+00 


a = 


ag{u)e 


iujO 


K 


c k e 


ikS 


k=—K 

K / +oo 


J2 c m 7, o(^) e 


,i(u+kLe)0 


k=—K \oj=— oo 

K / +oo 


) 


^2 Ck ( ^2 - kL e)( 

— -K \(jJ= — OO 

+oo / K \ 

E E c k a e {u - kLg ) 


iu)Q 


AujO 


uj=—oo \k=—K 


( 10 ) 


(ii) 


( 12 ) 


(13) 
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Then, (fTOl) becomes 


mm 

“£>,C/c 


+OO 

E 


K 


fe(u) - ^ c k a e (uj - kLg ) 


k=—K 

subject to: a£V(0). 


AujO 


2,0 


Then, using the well known Parsarval equality, the objective function is equal to the l 2 norm of the 
Fourier coefficients, which gives rise to the following equivalent optimization problem 

2 


nun 

ae,Ck 


+oo 

E 


UJ= — OQ 


K 


fe(u) - ^2 c k a e {u - kL g ) 


k=—K 


subject to: a £ V(9). 

Since a G V(6), using the dehnation of V(9), we have ag(cu) = 0, |w| > Lg/2. Then 


+oo 

E 


K 


feA) - ^2 c kae(u - kL e ) 


k=—K 


fg(u + jLg) - ^2 °Ae (w + (j - k)Lg) 


+oo LqI 2—1 K 

= E E 

j=—oou)=—Lg/ 2 k=—K 

K Lg/2-1 

= ^2 feA + jLg) - c k ag(uj) 

j=-K ej=-Lg/2 
Lei 2-1 

+ JeA+jLg) 

\j\>K, w=—Lg /2 

je z 

where we have used the fact that if k ^ j, ag(oj + (j — k)Lg) = 0 for any —Lg/2 < u ; < Lg/2 in 
obtaining the last equality. 

Using the above derivation, we have the following equivalent optimization problem, 

K 


min E E \fg(u + kLg ) - c k ag(uj)\ 2 . 




k=—K \uj\<Lq/2 


Denote 


fe,k( u ) = 


fg(uj), kLg < uj < (k + 1 )Lg, 
0, otherwise, 


and 


+00 


fg, k (0) = Fg 1 (fg, k (A)(0)= J2 MAe iuid . 
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Then, using the Parsaval equality one more time, we need only to solve the following equivalent 
problem: 

K 

min V \\fe,k(0) - c k a e(0)\\l, (14) 

k=—K 

where ag(6(t)) = a(t ) is the representation of a in the 0-space. Fortunately, after discretization, 
the above optimization problem can be solved by singular value decomposition (SVD). 

Suppose ag and fg y k is sampled over 6j = (j — 1 )/N, j = 1, ■ ■ ■ , N which is a uniform grid in 
the 0-space. Let 

f e,k = {fe,k(0 i), • • • , IoAOn))* , (15) 

F e = (Re(fe j0 ), • • • , Re^*-)), Im(f (9j i), • • • , Im(f e ,x)), 


and 

a 9 = (a e (0i), ■ ■ ■ ,0^(0^))*, (16) 

c = (Re(c 0 ), • • • , Re(c^), Im(ci), • • • , Im(c^)). 

Then, the discrete version of (1141) is 

min ||F* - a e -c\\ 2 F , (17) 

c6K 2if + 1 , 

a. e m N 

where || • \\f is the Forbenius norm of a matrix. It is well known that the above optimization 
problem can be solved by SVD. 

Suppose 

Ffl = U • S • V 

is the singular value decomposition of F#, diag(5) = (si, • • • , S 2 K+ 1 ) and si > S 2 > • • • , > 52^+1 > 
0. Then the solution of the optimization problem (1171) is 

a <9 = Ui, c = vi, (18) 

where ui is the first column of matrix U and Vi is the first row of matrix V. 

Summarizing above discussion, we get Algorithm [T] to compute the shape function with given 
phase function. 

At the end of this section, we give some remarks on how to compute the phase function 0. 
Notice that s has the following representation 

K K 

s(t ) = ^ Cfce* fc< = ^ (6fc cos (kt) + d & sin (kt )), 
k=-K k= 1 
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Algorithm 1 (Extraction of shape function) 

Input: Signal f = ■ ■ ■ ,f(tN )) is sampled over ti, l = 1, ■■■, N, phase functions 9, band 

width of shape function K. 

Output: Shape function s and corresponding envelope a(t). 

1: Interpolate the original signal f from ti,l = 1, • • • , N to a uniform grid 9j = j / N , j = 0, • • • , N— 
1 in 9 space. 

fg = Interpolate^^), 


In the computation, we use cubic spline to do interpolation. 
2: Compute the Fourier coefficients of / in 9 space 


N-l 

fo(v) = ft 

3=0 


3 e -iui9j 


N 


uj — — 



3: Compute fgk = 0, • • • , K , 

(fc+l/2)L e -l 

fe,k(0j)= E 

oj=(k—1/2) Lq 


4: Assemble the matrix according to (1151) . 

5: Apply SVD on F# to get the envelope ag(9j),j = 1, • • ■ , N and c*;, k = 0, • • • , K according to 

m and (usd. 

6: Compute the shape function and the envelope function in t space 

K 

s(ti ) = E c k f ' ktl - l = 

k=—K 

a(ti ) = Interpol&te(9j,ag(6j),6(ti)), 


where C-k = c* k is the complex conjugate of c^. The interpolation is also implemented by cubic 
spline. 





where bk = Re(cfc), dk = —Im(cfc). Without loss of generality, we assume that cq = 0. Otherwise, 
the constant part of s can be absorbed into r(t) in model (}7|) . Then the signal f(t) can be written 
as follows: 


f(t) = a(t)s(0(t)) + r(t) 

K 

= «(t)E (bk cos (k6(t)) + dk sm(k9(t))) + r(t) 
k =1 

= X] + d l) cos (k0(t) + <t> k ) + r(t ), 

k =l ' ' 

where 4>k = arctan ■ 

Then the signal /(f) can be seen as the signal composed by K IMFs. And these IMFs satisfy 
the scale separation property. Then, we could use the method developed in [12] to compute the 
phase function. 


4 The signal with varying shape function 


In some problems, the change of the shape function is more useful than the shape function itself. 
To detect the change of the shape function, we need some local method to extract “instantaneous” 
the shape function. The simplest idea is to cut the whole signal into small pieces and apply the 
method proposed in the previous section to get the shape function in each piece. Then we can get 
a series of shape functions which could show us how the shape function varies. 

Suppose we have a signal f = (/(ti), • • • , /(fjv )) which is sampled at t\, ■ ■ ■ ,tx, and the phase 
functions 8 = (8(t\),--- ,#(tjv)) is also given. For each t m , m = 1, • • ■ , N, we want to use the 
signal around t m to get a shape function. 

First, we extract a small piece of signals f m and the phase function 9 m around t m , the length 
of f m depends on the phase function 8, 


fm = frXT, 8 m = d T - 


where fr = ( f(tj))j&T and so is 8t- Here T = {1 < j < N : \9(tj) — 8(t m ) \ < /i7r} and 


XT = 



^1 + cos 



Here fx is a parameter to control the length of the segment. In this paper, we choose /i = 3, which 
means that for each point, we localize the signal within 3 periods to extract the shape function. 

Once we get the segments f m and 8 m , the shape function can be obtained by using the method 
in the previous section. For each t m , repeat this process, then we get a series of shape functions 
which could capture the change of the shape function. 
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The shape functions are a series of functions. Usually, it is not easy to distinguish which 
part is changing. It would be very helpful if we could find an index for each shape function 
and this index could reflect the main feature of the shape function. The concept of degree of 
nonlinearity was proposed by Huang (lecture in the IMA Hot Topic Workshop on Trend and 
Instantaneous Frequency, September 7-9, 2011, IMA). The main idea is that the signal comes from 
some physical process. The main feature of the shape function is the nonlinearity of the underling 
physical process. He defined an index to measure this nonlinearity. Later, we further developed this 
idea by assuming that the underling process is governed by a second order ODE with polynomial 
nonlinearity unj. Then we formulated an optimization problem to calculate the coefficients and 
the degree of nonlinearity of the ODE. Using these techniques, we could see clearly how the shape 
function changes and detect the time when significant change occurs. 


5 Numerical results 


In this section, we will present some numrical results to demonstrate the performance of our algo¬ 
rithm. 

Example 1 . The first example is a simple synthetic data which generated by the following formula: 


m 

a{t) 

/CO 


407rt + 2 cos( 67 rf), 

1 

2 + sin(27rt) ’ 

_a(t)_ 

1.1 + cos [9(t) + cos(20(t))] 


(19) 


From (1191) . the exact frequency is 0'(t) — 26'{t) sin(20(t)) which has two times more oscillations 
than the original signal f(t). If we relax the shape function to any periodic function not necessarily 
cosine, then the phase function 6 is much smoother than f(t) and the shape function has the form 
cos(f + cos(2f)). In Fig. [T] the shape function given by our method is shown. We can see, in this 
case, our method could capture the shape function very accurately. Next, we add Gaussian noise, 
0.3 X(t), to the clean signal to test the robustness of our method. X(t) is the standard Gaussian 
noise. The result is shown in Fig. [2j Even with noise, our method still recovers the shape function 
with reasonable accuracy. 


Example 2. The second example is the solution of the Duffing equation. This is an example 
to demonstrate the importance of the intra-wave frequency modulation in some complex dynamic 
system. 

The Duffing equation is a nonlinear ODE which has the following form: 

+ u + eu 1+UJ = 7 cos (f3t). (20) 
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Figure 1: Upper: the original data in Example 1; Bottom: The shape function obtained by our 
method (blue) and the exact shape function (red). 



Figure 2: Upper: The noised data f(t ) + 0.3X(f), where f(t ) is given in Example 1 and X(t ) is the 
white noise with standard derivative a 2 = 1; Bottom: The shape function obtained by our method 
(blue) and the exact shape function (red). 


11 





















































































Figure 3: Upper: the solution of duffing equation; Middle: the shape function s; Bottom: the 
Fourier coefficients of s. 

The parameters, e, 7 , w, that we use here to generate the solution in Fig. [3l are the same as those 
in the paper [7J, e = —1, 7 = 0.1, P = and oj = 2. The initial condition is u(0) = u'(0) = 1. 

In Fig. [3l we plot the shape function that we obtain from the solution of the Duffing equation. In 
[7], the example of the Duffing equation was used to demonstrate that the EMD method is capable 
of capturing the intra-wave frequency modulation. In that computation, the shape function is 
actually fixed to be the cosine function and the intra-wave oscillation is reflected in the instantaneous 
frequency. In our method, since the intra-wave oscillation is absorbed in the shape function, the 
instantaneous frequency is very smooth, but the shape function is not the simple cosine function 
any more, see Fig. [3j In this example, we can also express s(9) in terms of cos Ok, from which 
we can recover the instantaneous frequency with intra-wave modulation. More interestingly, from 
the Fourier coefficients of the shape function, we can see the Fourier coefficient is 0 when the 
wavenumber k = 2. This actually implies that u = 2 in the Duffing equation if we know that the 
signal comes from an ODE of the particular form given in (1201) . This phenomena may be very 
special, but it suggests that the some quantities, such as the deviation of the Fourier coefficients of 
the shape function, may reflect some important feature of the shape function. 

We also add Gaussian noise X(t) with variance a 2 = 1 to the original solution of the Duffing 
equation. Fig. []] shows the corresponding results. We can see that the shape function extracted 
from the noisy signal still keeps the main characteristics of the shape function extracted from the 
signal without noise. We can also clearly extract the degree of nonlinearity, w = 2, even with such 
large noise perturbation to the solution of the Duffing equation. 

Example 3: ECG data The last example is a piece of electrocardiogram (ECG) data. The 
length of the data used here is 16s. The bottom picture of Fig. [5]shows the shape function extracted 
from this set of ECG signal. We want to remark that it is challenging to extract the shape function 
from ECG data since it has sharp peaks in each period. This means that the shape function is 
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Figure 4: Upper: the solution of duffing equation with noise X(t)\ Middle: the shape function s; 
Bottom: the Fourier coefficients of s. 



Figure 5: Left: The original ECG data; Right: The shape function obtain by our method for the 
ECG data. 

not regular and needs many Fourier coefficients to well represent it. The shape function that we 
extracted seems to have all the characteristics of a typical ECG period. The interpretation of the 
significance of the shape function requires expertise in medicine and is beyond our expertise. 

6 Concluding Remarks 

In this paper, we present an effective and efficient method to extract the shape function from 
the signal with intra-wave frequency modulation by exploiting the intrinsic low rank structure of 
the data. The current method works only for those signal with one dominated shape function. 
Extracting shape functions for signals with multiple shape functions is much more involved and 
requires more efforts. How to define an effective index, such as the degree of nonlinearity, to reflect 
the main characteristic of the shape function is another interesting problem, a topic for our future 
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